//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : wvin_comp_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef WVIN_COMP_KERNELS_H
#define WVIN_COMP_KERNELS_H

namespace incomingwaves
{

template <class T>
__global__ void kernel_wvin_comp(T* zeta_in_X_loc,
                                 T* phi_in_X_loc,
                                 T* psi_in_X_loc,
                                 const int* wvin_X_count,
                                 const T* wvin_in_X_loc_ampl,
                                 const T* wvin_in_X_loc_phase,
                                 const T* wvin_in_X_loc_freq,
                                 const T* wvin_in_X_loc_facphi,
                                 const T* wvin_in_X_loc_facpsi,
                                 const int length,
                                 const T t,
                                 const T pslow,
                                 const int number)
{
    int tid = bx * Bx + tx;
    int idx = tid;

    T zeta = 0.0;
    T phi  = 0.0;
    T psi  = 0.0;

    // loop over waves (each thread does more work)
    //for (int wv = 1; wv <= number; ++wv) // fixed number for all threads
    for (int wv = 1; wv <= wvin_X_count[tid]; ++wv)
    {
        T phase  = wvin_in_X_loc_phase[idx] - wvin_in_X_loc_freq[idx] * t;
        T ampl   = wvin_in_X_loc_ampl [idx];
        T facphi = wvin_in_X_loc_facphi[idx];
        T facpsi = wvin_in_X_loc_facpsi[idx];

        T sinval;
        T cosval;

        // fastest alternative (look-up tables are not faster), reduces
        // computation time by 30% compared to sin() and cos() separately
        // however, still a rather expensive operation -> compute bound kernel
#ifdef USEDOUBLE
        sincos(phase, &sinval, &cosval);
        //__sincos(phase, &sinval, &cosval);
#else
        sincosf(phase, &sinval, &cosval); // slightly more accurate
        //__sincosf(phase, &sinval, &cosval); // slightly faster
#endif

        zeta += ampl * cosval;
        phi  += facphi * ampl * sinval;
        psi  += facpsi * ampl * sinval;

        idx  += length;
    }

    // write results to global memory
    zeta_in_X_loc[tid] = pslow * zeta;
     phi_in_X_loc[tid] = pslow * phi;
     psi_in_X_loc[tid] = pslow * psi;
}

} // namespace incomingwaves

#endif // WVIN_COMP_KERNELS_H

